#ifndef JETCLUSTERING_H
#define JETCLUSTERING_H
#include "utility.h"
#include "observables.h"
#include <stdio.h>
#include <iostream>
#include "mt2_bisect.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TF1.h"
#include "TUnuran.h"
#include "TUnuranContDist.h"

using namespace std;
using namespace mt2_bisect;

//*************************************************
extern "C" vector<vector<int> > jade(vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, vector<int> tobeclustered,
		float mcut) {
	//*************************************************

	vector<vector<int> > clusters;
	clusters.resize(0);
	/* Fill Clusters with a nested copy of tobeclustered */
	for (Int_t pcount = 0; pcount < tobeclustered.size(); ++pcount) {

		vector<int> temp;
		temp.resize(0);
		temp.push_back(tobeclustered[pcount]);
		clusters.push_back(temp);

	}

	bool goon;
	goon = true;

	while (goon) {

		int ntobeclulstered;
		ntobeclulstered = clusters.size();

		if (ntobeclulstered == 1) {
			goon = false;
		} else {

			vector<vector<vector<int> > > allkplets;
			vector<vector<int> > allpairs;

			findkplets(allkplets, 2, ntobeclulstered);

			allpairs = allkplets[0];

			//printmatrix(allpairs);

			//---->loop over *pairs*<------

			vector<int> minpair;
			minpair.resize(0);
			double minmassinv = 0;

			for (Int_t paircount = 0; paircount < allpairs.size(); ++paircount) {

				vector<vector<int> > sumofclusters;
				sumofclusters.resize(0);

				//		cout << "i vettori da mergiare" << endl;
				//		printvector(clusters[allpairs[paircount][0]-1]);
				//		printvector(clusters[allpairs[paircount][1]-1]);

				//sumofclusters = mergevector(clusters[allpairs[paircount][0] - 1], clusters[allpairs[paircount][1] - 1]);

				sumofclusters.push_back(clusters[allpairs[paircount][0] - 1]);
				sumofclusters.push_back(clusters[allpairs[paircount][1] - 1]);

				//		printvector(sumofclusters);
				double currentmass;

				currentmass
						= minvofn(sumofclusters, vec_partvector, partvector);
				//cout << "invmass " << currentmass << endl;

				if (1 / currentmass > minmassinv) {

					minmassinv = 1 / currentmass;

					minpair = allpairs[paircount];

				}

			}
			//			puts("****step*****");
			//			printmatrix(clusters);
			//			cout << "min invmas=" << 1.0 / minmassinv << "; minpair="; printvector(minpair); cout << endl;

			if (mcut > 1.0 / minmassinv) {

				if (minpair[0] > minpair[1]) {

					int temp;
					temp = minpair[0];
					minpair[0] = minpair[1];
					minpair[1] = temp;

				}

				vector<int> newclusters;
				int elem;
				elem = minpair[0] - 1;

				//printvector(clusters[ elem ]);

				newclusters = mergevector(clusters[minpair[0] - 1],
						clusters[minpair[1] - 1]);
				//				cout<< "il nuovo cluster " ; printvector(newclusters);
				//				cout<<endl;
				clusters = stripelement(minpair[1] - 1, clusters);
				clusters = stripelement(minpair[0] - 1, clusters);

				clusters.push_back(newclusters);

				//clusters e` stato aggiornato

				// printmatrix(clusters);

			} else {
				goon = false;
			}
		}

	}

	return clusters;
}

//*************************************************
extern "C" vector<vector<int> > kT(vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, vector<int> tobeclustered,
		float mcut) {
	//*************************************************

	vector<vector<int> > clusters;
	clusters.resize(0);
	/* Fill Clusters with a nested copy of tobeclustered */
	for (Int_t pcount = 0; pcount < tobeclustered.size(); ++pcount) {

		vector<int> temp;
		temp.resize(0);
		temp.push_back(tobeclustered[pcount]);
		clusters.push_back(temp);

	}

	bool goon;
	goon = true;

	while (goon) {

		int ntobeclulstered;
		ntobeclulstered = clusters.size();

		if (ntobeclulstered == 1) {
			goon = false;
		} else {

			vector<vector<vector<int> > > allkplets;
			vector<vector<int> > allpairs;

			/* ********************************************************************* */
			/* HERE I CAN READ THE FILE pairs_in_N.root TO GET allkplet[0] or allpairs*/
			/* ********************************************************************* */

			//			TFile *filepair;
			//			filepair = TFile::Open("pairs_in_N.root", "READ");

			//			allpairs = getthepairsfromfile(ntobeclulstered, filepair);
			allpairs = getthepairs(ntobeclulstered);

			//			filepair->Close("R");
			//			delete filepair;

			/* ********************************************************************* */
			/* AT THIS POINT I ALREADY HAVE allpairs*/
			/* ********************************************************************* */

			//printmatrix(allpairs);

			//---->loop over *pairs*<------

			vector<int> minpair;
			minpair.resize(0);
			double mininversekTmeasure = 0;

			for (Int_t paircount = 0; paircount < allpairs.size(); ++paircount) {

				vector<vector<int> > sumofclusters;
				sumofclusters.resize(0);

				//		cout << "i vettori da mergiare" << endl;
				//		printvector(clusters[allpairs[paircount][0]-1]);
				//		printvector(clusters[allpairs[paircount][1]-1]);

				//sumofclusters = mergevector(clusters[allpairs[paircount][0] - 1], clusters[allpairs[paircount][1] - 1]);

				sumofclusters.push_back(clusters[allpairs[paircount][0] - 1]);
				sumofclusters.push_back(clusters[allpairs[paircount][1] - 1]);

				//		printvector(sumofclusters);
				double currentkTmeasure;

				currentkTmeasure = kTmeasureofn(sumofclusters, vec_partvector,
						partvector);
				//minvofn(sumofclusters, vec_partvector, partvector);
				//cout << "invmass " << currentmass << endl;

				if (1 / currentkTmeasure > mininversekTmeasure) {

					mininversekTmeasure = 1 / currentkTmeasure;

					minpair = allpairs[paircount];

				}

			}
			//			puts("****step*****");
			//			printmatrix(clusters);
			//			cout << "min invmas=" << 1.0 / minmassinv << "; minpair="; printvector(minpair); cout << endl;

			if (mcut > 1.0 / mininversekTmeasure) {

				if (minpair[0] > minpair[1]) {

					int temp;
					temp = minpair[0];
					minpair[0] = minpair[1];
					minpair[1] = temp;

				}

				vector<int> newclusters;
				int elem;
				elem = minpair[0] - 1;

				//printvector(clusters[ elem ]);

				newclusters = mergevector(clusters[minpair[0] - 1],
						clusters[minpair[1] - 1]);
				//				cout<< "il nuovo cluster " ; printvector(newclusters);
				//				cout<<endl;
				clusters = stripelement(minpair[1] - 1, clusters);
				clusters = stripelement(minpair[0] - 1, clusters);

				clusters.push_back(newclusters);

				//clusters e` stato aggiornato

				// printmatrix(clusters);

			} else {
				goon = false;
			}
		}

	}

	return clusters;
}

//*************************************************
extern "C" vector<vector<int> > mcclusters(
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector, vector<int> tobeclustered,
		float mcut, int tipo) {
	//*************************************************

	//type=0,1,-1 select if the origintating partons are only the positive, the negative pids, any or both

	int type = tipo;

	// clusters will be the output
	vector<vector<int> > clusters;
	clusters.resize(0);

	vector<vector<int> > wipclusters;
	wipclusters.resize(0);

	//cout<<"/* COUNT THE MOTHERS BY A LOOP ON THE FULL EVENT*/"<<endl;
	int nmothers;
	nmothers = 0;

	for (int iparticle = 0; iparticle < partvector.size(); iparticle++) {
		if ((partvector[iparticle]->Status == -23)
				|| (partvector[iparticle]->Status == -24)) {
			nmothers++;
		}
	}

	//cout<<"no -23 or -24 found"<<endl;

	/* SCAN THE PARTICLE LIST ONE TIME PER EACH MOTHER AND ADD HER CHILDREN TO THE SAME CLUSTER */
	/* the Py8 LHE already know each final state parton from which mother was originated*/

	int steps = 1;
	if (type == 97) {
		steps = 2;
		type = 95;
	}

	// type 97 is imlpemendted by first doing the type 95 and then doing the type 96


	for (int cycle = 0; cycle < steps; cycle++) {

		for (int jmother = 0; jmother < nmothers; jmother++) {

			bool suitable;
			suitable = false;

			if (type == 94) {
				suitable = true;
			}

			if (type == 95 && partvector[jmother]->PID > 0) {
				suitable = true;
			}

			if (type == 96 && partvector[jmother]->PID < 0) {
				suitable = true;
			}

			if (suitable) {
				vector<int> thiscluster;
				thiscluster.resize(0);
				for (int nfinalpartons = 0; nfinalpartons
						< tobeclustered.size(); nfinalpartons++) {

					if (partvector[tobeclustered[nfinalpartons]]->Mother1
							== jmother + 1)
						thiscluster.push_back(tobeclustered[nfinalpartons]);

				}
				clusters.push_back(thiscluster);

			}
		}

		type++;
	}

	/*
	 * if we deal with pythia6 output no status -23 or -24 is present and charge of the jets have to be reconstucted summing all the charges
	 *
	 */

	/** for pythia6 output I assume the status is correcly set to "1 to 10" in the analysis.human file then all the parton are good **/

	/* find all the partciles with common mother
	 *     1) take the first particle mother index
	 *     2) add it to list of used mothers
	 *     3) take the next particle and read its mother index
	 *     4) check if it is the first time this mother is seen and in case add it to the list of used mothers
	 *     5) add this particle to the cluster corresponding to its mother
	 *
	 */

	cout<<"pythia jets"<<endl;

	if ((type == 98) || (type == 99) || (type == 100)) {
		vector<int> usedmothers;
		for (int nfinalpartons = 0; nfinalpartons < tobeclustered.size(); nfinalpartons++) {

			int thismother;
			thismother = partvector[tobeclustered[nfinalpartons]]->Mother1;

			if (not vectorcontains(usedmothers, thismother)) {

				usedmothers.push_back(thismother);
				/* correspondingly there is no cluster for this mother */
				vector<int> clusterofoneparticle;
				clusterofoneparticle.push_back(tobeclustered[nfinalpartons]);
				wipclusters.push_back(clusterofoneparticle);

			}

			/*
			 *  at this stage either the mother as been added to the list or was already there
			 * in any case we can look for it in the list and get its position
			 *
			 */
			int clusterindex;
			clusterindex
					= vectorcontainsatposition(usedmothers, thismother).second;
			wipclusters[clusterindex].push_back(tobeclustered[nfinalpartons]);

		}

		/* now clusters contains the clusters but it is not known which is the positive and which is negative charge */

		/*
		 * 0) check if a specific charge is required
		 * 1) determine the charge of each cluster
		 * 2) drop the unwanted component
		 */
		if (type == 100) {
			clusters = wipclusters;
		}

		if (type == 98) {

			cout<<"pythia jets type 98"<<endl;

			for (int jcluster = 0; jcluster < wipclusters.size(); jcluster++) {
				float charge;
				vector<vector<int> > matrixtoadd;
				matrixtoadd.push_back(wipclusters[jcluster]);

				charge
						= totalchargeofn(matrixtoadd, vec_partvector,
								partvector);
				//find charge of the cluster
				if (charge < 0.) {
					clusters.push_back(wipclusters[jcluster]);
				}

			}
		}

		if (type == 99) {
			//find charge of the cluster
		}

	}

	return clusters;

}

#endif
